home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / FIT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  2KB  |  64 lines
  1. PROCEDURE fit(x,y: glndata; ndata: integer; sig: glndata; mwt: integer;
  2.          VAR a,b,siga,sigb,chi2,q: real);
  3. (* Programs using routine FIT must define the type
  4. TYPE
  5.    glndata = ARRAY [1..ndata] OF real;
  6. in the main routine.   *)
  7. VAR
  8.    i: integer;
  9.    wt,t,sy,sxoss,sx,st2,ss,sigdat: real;
  10. BEGIN
  11.    sx := 0.0;
  12.    sy := 0.0;
  13.    st2 := 0.0;
  14.    b := 0.0;
  15.    IF (mwt <> 0)THEN BEGIN
  16.       ss := 0.0;
  17.       FOR i := 1 TO ndata DO BEGIN
  18.          wt := 1.0/sqr(sig[i]);
  19.          ss := ss+wt;
  20.          sx := sx+x[i]*wt;
  21.          sy := sy+y[i]*wt
  22.       END
  23.    END ELSE BEGIN
  24.       FOR i := 1 TO ndata DO BEGIN
  25.          sx := sx+x[i];
  26.          sy := sy+y[i]
  27.       END;
  28.       ss := ndata
  29.    END;
  30.    sxoss := sx/ss;
  31.    IF (mwt <> 0)THEN BEGIN
  32.       FOR i := 1 TO ndata DO BEGIN
  33.          t := (x[i]-sxoss)/sig[i];
  34.          st2 := st2+t*t;
  35.          b := b+t*y[i]/sig[i]
  36.       END
  37.    END ELSE BEGIN
  38.       FOR i := 1 TO ndata DO BEGIN
  39.          t := x[i]-sxoss;
  40.          st2 := st2+t*t;
  41.          b := b+t*y[i]
  42.       END
  43.    END;
  44.    b := b/st2;
  45.    a := (sy-sx*b)/ss;
  46.    siga := sqrt((1.0+sx*sx/(ss*st2))/ss);
  47.    sigb := sqrt(1.0/st2);
  48.    chi2 := 0.0;
  49.    IF (mwt = 0)THEN BEGIN
  50.       FOR i := 1 TO ndata DO BEGIN
  51.          chi2 := chi2+sqr(y[i]-a-b*x[i])
  52.       END;
  53.       q := 1.0;
  54.       sigdat := sqrt(chi2/(ndata-2));
  55.       siga := siga*sigdat;
  56.       sigb := sigb*sigdat
  57.    END ELSE BEGIN
  58.       FOR i := 1 TO ndata DO BEGIN
  59.          chi2 := chi2+sqr((y[i]-a-b*x[i])/sig[i])
  60.       END;
  61.       q := gammq(0.5*(ndata-2),0.5*chi2)
  62.    END;
  63. END;
  64.